close all;clear all;clc;
% Parameters
A = 1.1;
g = 2;
lambda = 0.2;
eta = 0.6;
gprime = 1;

% Functions
tauf = @(g, gprime, lambda, eta, A, k)...
    1/(g-gprime)*(log(A)-log((exp((g-gprime)*k*eta)-1)/((g-gprime)*eta)...
    + (1-k)*((g/lambda)*exp((g-lambda-gprime)*k*eta)...
    - (g-lambda)/lambda*exp((g-gprime)*k*eta))));

alphaf = @(g, gprime, eta, A, k, tau)...
    (A*exp(gprime*(tau+k*eta))...
    - exp(g*tau)*(exp(g*k*eta)-exp(gprime*k*eta))/((g-gprime)*eta))/(1-k);

hf = @(g, lambda, eta, k) lambda*exp(lambda*k*eta)/(exp(lambda*k*eta)-1);

kcf = @(g, gprime, lambda, eta) fsolve(@(k) (1-k)*g*eta + g/lambda...
    - (g-lambda)/lambda*exp(lambda*k*eta) - exp(-(g-lambda-gprime)*k*eta), 0);

% Calculations
kx = 0:0.001:1/((g-gprime)*eta)*log(A*(g-gprime)*eta+1);

for i = 1:length(kx)
    taux(i) = max(tauf(g, gprime, lambda, eta, A, kx(i)),0);
    if taux(i)>0
        k(i) = kx(i);
        tau(i) = taux(i);
        alpha(i) = alphaf(g, gprime, eta, A, k(i), tau(i));
        F(i) = (exp(lambda*eta) - exp(lambda*k(i)*eta))/(exp(lambda*eta) - 1);
        h(i) = hf(g, lambda, eta, k(i));
        
        tauPlusKeta(i) = tau(i) + k(i)*eta;
        decomp1_h(i) = 1/g*log(h(i)/(h(i)-g));
        decomp2_keta(i) = k(i)*eta;
        decomp3_alpha(i) = 1/g*log(alpha(i)*exp(-g*k(i)*eta));
      
        diff(i) = tauPlusKeta(i) - (decomp1_h(i)+decomp2_keta(i)+decomp3_alpha(i));
        decomp1_h_end = decomp1_h(i);
        decomp3_alpha_end = decomp3_alpha(i);
    end
    
    if taux(i) == 0
        k(i) = kx(i);
        tau(i) = taux(i);
        
        tauPlusKeta(i) = tau(i) + k(i)*eta;
        decomp1_h(i) = nan;
        decomp2_keta(i) = k(i)*eta;
        decomp3_alpha(i) = nan;
    end
end

knc = min(1/(lambda*eta)*log(g/(g-lambda)),...
    1/((g-gprime)*eta)*log(A*(g-gprime)*eta+1));
kc = kcf(g, gprime, lambda, eta);
ykc = tauf(g, gprime, lambda, eta, A, kc) + kc*eta;

% Points for figures
for j = 1:length(k)
    if (j==1)|(j==length(k))|(mod(j,50)==0)
        k_figure(j) = k(j);
        tauPlusKeta_figure(j) = tauPlusKeta(j);
        decomp3_alpha_figure(j) = decomp3_alpha(j);
        decomp1_h_figure(j) = decomp1_h(j);
        decomp2_keta_figure(j) = decomp2_keta(j);
    else
        k_figure(j) = NaN;
        tauPlusKeta_figure(j) = NaN;
        decomp3_alpha_figure(j) = NaN;
        decomp1_h_figure(j) = NaN;
        decomp2_keta_figure(j) = NaN;
    end
end

% Figures
h = figure;
set(h, 'defaultLegendInterpreter','latex',...
    'defaultAxesTickLabelInterpreter','latex');
a = axes;

plot(k,tauPlusKeta,'k','linewidth',2);
hold on;
plot(k,decomp3_alpha_figure,'linewidth',2,...
    'Marker','*','LineStyle', '--','color', 'k', 'MarkerSize',8);
plot(k,decomp1_h_figure,'linewidth',2,...
    'Marker','o','LineStyle', '-.','color', 'k', 'MarkerSize',8);
plot(k,decomp2_keta_figure,'linewidth',2,...
    'Marker','square','LineStyle', ':','color', 'k', 'MarkerSize',8);
yl = ylim();
plot([kc kc], [yl(1) ykc], 'k--', 'linewidth', 1);

plot(k,decomp3_alpha,'LineStyle', '--','color', 'k','linewidth',1);
plot(k,decomp1_h,'LineStyle', '-.','color', 'k','linewidth',1);
plot(k,decomp2_keta,'LineStyle', ':','color', 'k','linewidth',1);

set(a, 'xlim', [0 knc], 'XTick',[0 kc knc],...
    'XTicklabel',{'$0$' '$k_c$' '$k_{nc}$'},...
    'ylim', [decomp3_alpha_end-0.05, decomp1_h_end+0.05],...
    'Fontsize', 16);
xlabel('Bankruptcy Threshold ($k$)', 'interpreter','latex');
legend1 = legend('$\tau^{*}(k)+k\eta$',...
    '$\frac{1}{g}\log\left[\alpha(\tau^{*},k) e^{-g k \eta}\right]$',...
    '$\frac{1}{g}\log\left(\frac{h(k)}{h(k)-g}\right)$',...
    '$k\eta$');
set(legend1,'Location','southwest');

%saveas(h, 'tauplusketa_decomposition_k.eps');



